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ABSTRACT 

We describe a data reduction pipeline for VLBI astrometric observations of 
pulsars, implemented using the ParselTongue AIPS interface. The pipeline per- 
forms calibration (including ionosphere modeling), phase referencing with proper 
accounting of reference source structure, amplitude corrections for pulsar scintil- 
lation, and position fitting to yield the position, proper motion and parallax. The 
optimal data weighting scheme to minimize the total error budget of a parallax 
fit, and how this scheme varies with pulsar parameters such as flux density, is also 
investigated. The robustness of the techniques employed are demonstrated with 
the presentation of the first results from a two year astrometry program using 
the Australian Long Baseline Array (LBA). The parallax of PSR J1559-4438 is 
determined to be vr = 0.384 ±0.081 mas (lo"), resulting in a distance estimate of 
2600 pc which is consistent with earlier DM and HI absorption estimates. 

Subject headings: Techniques: interferometric — astrometry — pulsars: general 
— pulsars: individual(J1559-4438) 



1. Introduction 

Accurate, model independent distances to pulsars are, like many astronomical distance 
measurements, highly prized but difficult to obtain. The dispersion measure (DM) of a 
pulsar indicates the integrated column density of free electrons between the observer and the 
pulsar, and can be used in conjunction with Galactic electron distribution models to estimate 



^Centre for Astrophysics and Supercomputing, Swinburne University of Technology, Mail H39, P.O. Box 
218, Hawthorn, VIC 3122, Austraha 

^co-supervised through the Australia Telescope National Facility, P.O. Box 76, Epping, NSW 1710, Aus- 
tralia 

■^Department of Imaging and Applied Physics, Curtin University of Technology, Bentley, WA, Australia 
^National Radio Astronomy Observatory, Charlottesville, VA, USA 



- 2 - 



pulsar distances, but this approach suffers from considerable uncertainty when the electron 
distribution model is poorly understood, such as at high Galactic latitudes. Timing residuals 
for millisecond pu lsars (MSPs) can be used to determine parallax and proper motion (e.g. 



Hotan et al.ll2006l ). but only for the very limited subset of pulsars with extremely accurate 
timing solutions. VLBI astrometry offers a means to directly measure the parallaxes and 
proper motions of nearby pulsars. 

By obtaining an independent pulsar distance estimate, a pulsar DM can be used to 
determine the average electron density along the line of sight; an ensemble of such mea- 
surements can be used to improve electron distribution models. Whilst several recent large 
puls ar parallax programs h ave significantly increased the number of known VLBI parallaxes 
(e.g. IChatterjee et al.ll2004l ). to date very few VLBI parallaxes have been obtained for south- 
ern pulsars, due to the bias of VLBI facilities towards the Northern Hemisphere. As such, 
models of Galactic electron distributions are more uncertain at far southern declinations. 

Furthermore, binary pulsars offer the opportunity to make exceedingly precise tests of 
gravitational theories, since relativistic effects contribute to the observed rate of change of 
binary period (Ph), which c an be measured very precisely in systems with a stable pulsar. 



However, kinematic effects (IShklovskiil Il970l ) also contribute to Ph, and can only be accu- 
rately subtracted to obtain P^ due to General Relativity (GR) if the pulsar distance and 
transverse velocity are accurately known. The effect scales with the square of proper mo- 
tion, and hence is typically largest for nearby pulsars, where VLBI measurements of parallax 
are most feasible. Since the uncertainty of a DM-based distance estimate is large for any 
individual pulsar, a reliable estimate of the error on derived GR quantities requires a direct 
distance determination. Similarly, luminosities for individual pulsars based on DM-inferred 
distances are questionable, so deriving accurate an accurate luminosity for any individual 
pulsar generally requires an independent confirmation of distance. 

In this paper, we describe our target selection policy in §2], and the observations un- 
dertaken in §31 The data reduction pipeline is described in SjH and we present the results 
of PSR J1559-4438 in f|6l We analyse the optimum visibility weighting scheme to use for 
VLBI astrometry in §3 and the magnitudes of different sources of systematic errors are 
investigated in §H1 Our conclusions are presented in ^ 



2. Target selection 



Our observational program encompassed 8 pulsars, listed in Table [H Previous Southern 
Hemisphere VLBI pulsar astrometry programs (jPodson et al.ll2003uLeggell2002l ; lBailes et al. 
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19901 ) have yielded only two published pulsar parallaxes, so one of our primary motivations 
was to increase the number of southern VLBI parallaxes and hence improve models of Galac- 
tic electron distributions at southern declinations. Additionally, as this is the first such large 
scale southern parallax study, we chose to target pulsars of varying brightness and predicted 
distance to determine the types of targets that would be feasible for future southern VLBI 
studies. 

In this paper we present results for the bright pulsar J1559-4438. As an isolated pulsar 
with a characteristic age of 4 Myr and DM distance of 2.35 kpc, J1559-438 is unremarkable, 
but as it provides the highest signal to noise detections of our target pulsars, and suffers 
most heavily from scintillation, it provides the clearest example of the techniques which we 
will apply to the data reduction of the remaining pulsars. PSR J2048-1616 was included as a 
second "technique check" source and those results will be presented in a future pubhcation. 

The remaining targets can be broadly divided into two scientific categories: binary pul- 
sars used for tests of GR and gravitational wave detection (J0437-4715, J0737-3039A/B 
and J2145-0750), and pulsars with unusual DM-based luminosity in the radio (low lumi- 
nosity pulsars J0108-1431 and J2144-3933) or x-ray (J0630-3834, whose x-ray luminosity 
is anomalously high). Each group will be the subject of a forthcoming analysis. 



3. Observations 

Eight observational epochs were spread over a two year period between May 2006 and 
February 2008. Epochs were typically 24 hours in duration, and observed a subset of the 
8 pulsars, depending on which were closest to parallax extrema. The observing frequency 
was centred on 1400 MHz for the first observation, and 1650 MHz for the remaining 7 obser- 
vations. PSR J0437-4715 was observed separately, with four epochs of 12 hours duration, 
centred on 8400 MHz. Observations at this higher frequency were made possible by the high 
flux and narrow pulse profile of PSR J0437-4715 . Dual circular polarization was used at all 
frequencies. 

The Australian Long Baseline Array (LBA) consists of six antennas - the Australia Tele- 
scope National Facility (ATNF) telescopes in New South Wales (Parkes, Australia Telescope 
Compact Array [ATCA], Mopra); the University of Tasmania telescopes at Hobart, Tasma- 
nia and Ceduna, South Australia; and the NASA DSN facility at Tidbinbilla, Australian 
Capital Territory. The Parkes, phased ATCA, Mopra, and Hobart telescopes participated 
in all experiments, but a Tidbinbilla antenna (70m or 34m) was used only when available. 



-4- 



and Ceduna participated in observations of J0437-4715 onl}0. The maximum baseline length 
with Ceduna is 1700 km, and without Ceduna is 1400 km. Representative uv coverage at 
1650 MHz and 8400 MHz is shown below in Figured 

All observations used the recently introduced LBADR disk-based recording system 
(Phillips et al., in preparation). At the three ATNF observatories, the presence of two 
Data Acquisition System (DAS) units allowed a recording rate of 512 Mbps (8 x 16 MHz 
bands, Nyquist sampled at 2 bits), while the non-ATNF stations recorded at 256 Mbps (4 x 
16 MHz bands). For epochs where dual-polarization feeds were available at all antennas, two 
frequency bands were dropped at the non-ATNF stations, but as some of the NASA DSN 
feeds are single polarization only, 1650 MHz epochs featuring the 70m NASA antenna and 
8400 MHz epochs featuring the 34m NASA antenna instead retained a single polarization of 
all frequency bands. 

A phase reference cycle time of six minutes, apportioned equally to target and calibrator, 
was used for all observations. As the LBA consists of disparate antennas ranging up to 70m 
in diameter (and, when phased, the ATCA has an equivalent diameter of hundreds of metres 
for the purposes of calculating field of view) , it was not possible to utilise in- beam calibrators 
for any sou rces, unlike recent pulsar astrometric programs using the VLBA Jchatteri ee et al.l 



2001 



2005f ). 



The data w ere correlated using matched filtering oripulse profiles with the DiFX soft- 



ware correlator ( Deller et al.l 120071 ). producing RPFITqj format visibility data. Matched 
pulse profile filters (a more advanced pulsar 'gate') allow the maximum recovery of signal to 
noise when observing pulsars, by dividing the pulse into bins and appropriately weighting 
each bin by the expected signal strength before summation. Table [T] shows the predicted 
gain due to pulse profile filtering f or each target source. P ulsar ephemerides were obtained 



using the ATNF Pulsar Catalogue (IManchester et al.ll2005l ). with the exception of the double 



pulsar J0737-3039, which was periodically updated with the latest published ephemeris. Two 
second integrations and 64 spectral points per 16 MHz band were used for all observations. 



^Ceduna does not possess a 1600 MHz receiver, so could only participate in the higher frequency experi- 
ments 

^http://www. atnf.csiro.au/computing/software/rpfits. html 
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Data reduction 



Dat a reduction was perf ormed principally in AlP qij, using the py thon interface Parsel- 
Tongue (IKettenis et al.ll2006l ). The DIFMAP package (IShepherdlll997l ) was used for imaging 
and self calibration. The data reduction was implemented as a pipeline, with user interaction 
for imaging, editing of solution tables, and visibility flagging. All scripts used in the data re- 
duction process are available at http : //astronomy, swin.edu.au/~adeller/ software/scripts/. 
The individual stages of the pipeline are described below. 



4.1. Amplitude and weight calibration and flagging 

Amplitude calibration using the measured telescope system temperatures was carried 
out using the AIPS tasks APCAL and ANTAB. Flagging based on predicted or (when avail- 
able) logged telescope off-source times due to slewing or failures was applied using UVFLG, 
while the first and last 10 seconds of every scan was excised with the task QUACK. For the 
ATCA, the tied array infrastructure required flagging the first 3 correlator integrations (30 
seconds) of each scan with QUACK. The weight of each visibility point, which is set to a 
constant value when the RPFITS format data is loaded into AIPS, was initially scaled by 
the predicted baseline sensitivity using a ParselTongue script. The effect of data weighting 
is investigated further in §6.31 



4.2. Geometric model and ionospheric corrections 

At the low frequencies which are generally used for pulsar astrometry, io nospheric vari 



ations usually make the dominant contribution to systematic error (see e.g. iBrisken et al. 



20021 ) . Using ionospheric models based on Global Position System (GPS) data provided 



by the NASA Jet Propulsion Laboratorjfl, we corrected for phase variations due to the 
ionosphere using the task TECOR. This observing program roughly coincided with the solar 
minimum of 2006, and consequently the ionospheric variations were generally at a minimum. 

Total Electron Content (TEC) models suffer at southern declinations due to the rela- 
tively sparse distribution of GPS receivers at southern latitudes. Consequently, the derived 



•^httpi/Zwww. aoc.nrao.edu/aips 

''available from the Crustal Dynamics Data Information Systems (CDDIS) archive: 
ftp : / / cddis . gsfc . nasa.gov / pub / gps / products/ ionex / 
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ionospheric corrections are much less rehable for the LBA than for similar Northern Hemi- 
sphere instruments. The dispersive delay corrections generated by TECOR were inspected 
for each epoch, and the effectiveness of different TEC maps is investigated in §6.21 

While the observational program was underway, considerably more accurate station po- 
sitions were derive d for several LBA antennas us ing archival geodetic observations and the 



OCCAM software (ITitov. Tesmer fc Boehmll2004[ l. a dedicated 22 GHz LBA geodetic exper 



iment (Petrov et al., in preparation) and GPS measurements. Additionally, more accurate 
positions for some calibrators were published in the 5th VLBA Calibrator Survey (VCS5; 



Kovalev et al.ll2007l ). The visibilities were corrected to account for the revised positions using 
an AIPS SN table generated with the Wizardry feature of ParselTongue, which stored the 
difference between the initial and corrected geometric models. 

For some pulsars, their proper motions (> 100 mas yr~^) caused significant position 
shifts over the course of a 24 hour observation, comparable in some cases to the epoch 
positional accuracy. As the geometric model generation used in DiFX at the time of these 
observations could not account for proper motion, the visibility phases and uvw values 
were corrected in AIPS using a SN table generated by ParselTongue, interpolating between 
predicted postions for the pulsar at the start and end of the experiment. 



4.3. Fringe-fitting and amplitude calibration refinement 

Fringe fitting was performed using the AIPS task FRING, using a point source model, on 
the phase reference calibrator data for each target pulsar. Subsequently, a single structural 
model was produced for each calibrator using the combined datasets from all epochs. Each 
source was modeled using a dominant component (delta or narrow Gaussian) fixed at the 
phase center, and 0-2 secondary components which were allowed to vary in position. The 
flux of all components was allowed to vary. Typically, the variation between epochs of 
secondary component (s) flux were < 1% of peak image flux. The solutions were applied 
to each calibrator and the data averaged in frequency, exported to disk and loaded into 
DIFMAP. The calibrator model was loaded and several iterations of self-calibration and 
modelfitting performed. The difmap 'modelfit' command uses the Levenberg-Marquardt 
least-squares minimisation algorithm to fit the free model parameters to the visibility points, 
incorporating the visibility weights. The self calibration corrections were then written to disk 
as an AIPS SN table using the 'cordump' patch to difmap|f|. Additionally, data points flagged 
in DIFMAP were collated in a Wizardry script and converted into a flag flle suitable for the 



^http://astrononiy.swin. edu.au/~elenc/difniap-patches/ 
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AIPS task UVFLG. 

The amplitude corrections generated in this manner allowed compensation for imper- 
fectly measured system temperature values, and were also applied to visibility weights. For 
bands which could not be self-calibrated, due to only being observed by the three ATNF an- 
tennas, a correction was estimated based on the nearest band with self calibration solutions 
of the same polarization. The self calibration solutions were loaded into AIPS using TBIN, 
and applied to the target pulsars using CLCAL. 

The use of bandpass calibration was investigated but found to make insignificant dif- 
ference to the fitted target position. The LBA Data Acquisition System (DAS) utilitizes 
digital filtering and typical bandpass phase ripple was < 2°. When averaging in frequency, 
the lowest and highest 10% of the band was discarded and the central 80% of the band 
averaged with uniform weight assigned to each channel. 



4.4. Pulsar scintillation correction 



Nearby pulsars can suffer dramatic, and rapid, variations in visibility amplitude due to 
diffractive scintillation. The size of the scintillation pattern is typically much larger than 
the size of the Earth, and so the amplitude variations are essentially independent of baseline 
length. Maximal amplitude fluctuations (where t he rms is equ al to the mean flux) are seen 
for pulsars in the strong scattering regime (see e.g. lWalkeiill998l), which can be predi cted from 
Galactic electron distribution models. The NE2001 model (jCordes &: Lazid l2002l ) predicts 
that strong scintillation should be observed for J0630-2834, J1559-4438, J2048-1616, J2144- 
3933 and J2145-0750. As an example. Figure [2^ shows the variation of visibility amplitude 
with time for J1559-4438 on Tidbinbilla baselines over a 2 hour period. Observed scintillation 
parameters for target pulsars are shown in Table [2l Assuming the material responsible for 
the scintillation is turbulent, with a Kolmogorov distribution, the scintillation time Tscint and 
scintillation bandwidth Ssdnt can be scaled to the frequencies u sed in these observations with 



the relations Tgcint oc ly^''^ and -Bscint oc z/^'"^ (ICordes et al.lll986l ). 



Left uncorrected, the variation in visibility amplitude with time scatters power through- 
out the image plane, as shown in Figure [Jt, which shows the residuals for J1559-4438 after 
fitting a single point source to the visibilities shown in Figure [2^. To overcome this, a 
ParselTongue script was written to produce an AIPS SN table which would flatten visibihty 
amplitudes by averaging data over all sensitive baselines to 1/4 of the scintillation time, nor- 
malizing, and taking the square root to obtain an antenna based correction. The visibility 
weights were then scaled by the inverse of the square of the correction, upweighting points 
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when the amphtude was high and downweighting points of low significance. The effect of 
these corrections is shown in Figure [2]d and d, which show the visibihty amphtudes and im- 
age residuals respectively. In this example, the signal to noise ratio (SNR) of the detection 
is improved by a factor of two when the visibility amplitudes are corrected for scintillation, 
which ultimately reduces the error on the position determination by a corresponding factor. 



5. Positional determination and parallax fitting 

The calibrated visibilities for each pulsar were averaged in frequency, written to disk 
and loaded into DIFMAP. A single delta function model component was initialized at the 
peak of the dirty image, and the modelfit command used to obtain the best fit for the pulsar 
visibility data. Each observing band, as well as the combined dataset, was then imaged 
separately using uniform weighting (2 pixels, weights raised to the power —1) producing 
8 images which were saved and read back into AIPS using the task FITLD. Variations to 
the weighting scheme are discussed in §6.31 The AIPS task JMFIT was used to determine 
the pulsar position and formal (SNR-based) errors in the image plane. Systematic offsets of 
several (3-6) mas were observed at all epochs between the 4 bands that were only contributed 
to by the three ATNF antennas, and the 4 bands common to all antennas. The magnitude 
and direction of these offsets varied between epochs, and so the ATNF-only bands (which 
additionally had formal errors of approximately five times the other bands, primarily due to 
the shorter baselines) were dropped. We suspect that the systematic offsets were caused by 
the lack of an accurate amplitude self-calibration solution for these bands. 

For each pulsar, the best fit values of J2000 position (RA and declination), proper motion 
(RA and declination) and parallax were initially determined by iteratively minimising the 
error function calculated by summing the value of predicted minus actual position over all 
measurements, weighted by the individ u al me asurement errors. The iterative minimisation 



code used is described in iBrisken et al.l ( 120021 ). A reference time for the proper motion was 



chosen to be 31 Dec 2006 (MJD 54100) to minimize proper motion uncertainty contibutions 
to be position uncertainty. 

This approach yields error estimates for the 5 fitted parameters which are almost cer- 
tainly an underestimate, for two reasons: 

1. There are systematic errors, varying from band to band within an epoch (intra-epoch 
errors) such as the residual unmodeled differential ionosphere, bandpass effects etc, 
and varying from epoch to epoch (inter-epoch errors) caused by effects dependent 
on observing time, such as seasonal or diurnal ionospheric variations, and variations 



in refractive scintillation image wander (IRickettl Il990l ). These should increase the 
error on each individual measurement, but estimating the magnitude of the systematic 
component for each individual measurement is poorly constrained. 



2. Each measurement is assumed to be completely independent, whereas as noted above 
we expect correlated errors between measurements from the same epoch. In essence, 
this approach overestimates the number of independent measurements, lowering the 
reduced chi-squared and implying a better fit than the actual result. 

It is possible to make an estimate of the magnitude of the intra-epoch errors by comparing the 
scatter in fitted positions from the individual bands. Forming a weighted centroid position for 
each epoch utilising all measurements from that epoch allows an estimation of the likelihood 
that the measured points are consistent with that centroid, through the calculation of a 
reduced chi-squared value. If the reduced chi-squared value exceeds unity, the presence of 
unmodeled systematic errors can be inferred. 

Since we have no a priori knowledge of the systematic error distribution, we have cho- 
sen to allocate an equal systematic error to each measurement, and add in quadrature to 
the original measurement error. The errors in right ascension and declination are treated 
separately. This is necessarily an iterative procedure, since the weighted centroid will be 
altered by the addition of these systematic error estimates. In effect, this assumes a zero 
mean, gaussian distribution for the systematic errors. Although this is unlikely to be the 
true distribution, it is the most reasonable assumption available, and certainly more correct 
than assuming no systematic errors at all. The intra-epoch systematic error estimate for an 
epoch is obtained when the reduced chi-squared reaches unity. 

Once a single position measurement and error has been calculated for each epoch, the fit 
to position, parallax and proper motion can be re-calculated, and the reduced chi-squared 
of the fit inspected again. Since the number of measurements were reduced, the addition 
of systematic errors did not always lower the reduced chi-squared of the fit. If reduced 
chi-squared remained significantly greater than unity, we concluded that significant inter- 
epoch systematic errors remained. As with intra-epoch errors, the true error distribution is 
unknown, and so the inter-epoch error was apportioned equally between epochs. The errors 
in right ascension and declination were treated separately, and iterated until a reduced chi- 
squared of unity was obtained. 

Thus, the final astrometric dataset for each pulsar consisted of a single position mea- 
surement for each epoch, with a total error equal to the weighted sum of the individual 
band formal errors, added in quadrature to the estimates of intra-epoch and inter-epoch 
systematic errors. 
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For each pulsar, the robustness of error estimation after the inclusion of estimated sys- 
tematic contributions was checked by implementing a bootstrap technique. Bootstrapping, 
which involves repeated trials on samples selected with replacement from the population 
of measured position points, is a statistical technique allowing the estimation of parameter 
without complete knowledge of the underlying distribution JPressll2002h . In this in- 



errors 

stance, the original single band position measurements were taken as the population, and 
N samples were drawn, where N was the original number of measurements. Each bootstrap 
consisted of 10,000 such trials, and the parameter errors estimated from the variance of the 
resultant distributions. A minimum of 3 different epochs needed to be included to ensure a 
fit was possible - on the rare occasion that a trial did not satisfy this requirement, it was 
re-drawn. 

Thus, three sets of fitted parameters and errors were obtained for each pulsar - a 
"naive" result using the single-band positions, a bootstrap result, and a more conservative 
estimate which attempts to account for the impact of systematic errors (the "inclusive" fit). 
In general, the estimated magnitude of errors on fitted parameters increased through these 
three different schemes. Typically, the ratio in the errors on the inclusive fit to those on 
the naive fit ranged from 0.95 to 1.90. We feel that the final error values obtained from the 
inclusive fit are the most accurate estimation possible, and are inherently conservative. All 
quoted errors are la unless otherwise stated. 



6. Results for PSR J1559-4438 

6.1. Initial results 

Using the techniques described above, we obtained initial results for J1559-4438 (shown 
in Table [3]). The motion of J1559-4438 in right ascension and declination is shown in 
Figure El along with the fitted path according to the systematic-error weighted fit. The fit 
is clearly unsatisfactory, since it predicts a negative parallax (though consistent with zero). 
The final column of Table [3] shows that systematic errors far exceed the nominal single-epoch 
positional accuracies, and inspection of Figure [3] shows that the first epoch (MJD 53870) is 
markedly discrepant with the remaining epochs. 

As we show below, fine-tuning of the data reduction is required in order to obtain 
optimal results from each pulsar, in particular with regard to the details of the ionospheric 
corrections and the visibility weighting schemes using in imaging. We describe the steps 
taken for J1559-4438 in gOand S 



6.2. Ionospheric correction 



The obvious source of the large systematic errors present in the initial fit shown in 
§6.1l is the varying ionospheric correction between epochs. Accordingly, as a first check, the 
position fits for each epoch were recalculated after subtracting the differential ionospheric 
correction between PSR J1559-4438 and its phase reference - in effect, removing the applied 
ionospheric correction and leaving the data uncorrected for ionospheric effects. This was 
implemented using a ParselTongue script which made a two-point interpolation between 
adjacent calibrator scans to calculate the differential correction to the target (which was not 
absorbed into the fringe-fit), which was stored in a CL table and subtracted using the AIPS 
task SPLAT. The applied values were saved for later analysis. 

Typically, we would expect the largest ionospheric corrections when the angular dis- 
placement of the pulsar from the sun is small, since angular displacement and solar activity 
are the largest influence on TEC. The angular displacement at each epoch between the orig- 
inal fitted position and the position obtained when ionospheric correction was removed is 
presented in Table HJ and plotted against angular separation of the pulsar from the sun at 
the time of observation in Figure |H The revised astrometric fit obtained without ionospheric 
correction is plotted in Figure [51 

It is immediately apparent from Figure H] that the first epoch (MJD 53870) is discrepant 
in that the position shift due to ionospheric correction is unusually large, given the large 
angular separation from the Sun. As shown in Figure [5], this epoch becomes more consistent 
with the fit when the ionospheric correction is removed, but the third epoch (MJD 54057) 
becomes much more inconsistent. This is unsurprising, however, since this epoch had the 
smallest pulsar-Sun separation and the largest ionospheric corrections. 

To investigate whether the chosen ionospheric map was at fault, the first epoch was re- 
reduced with all available maps from the NASA CDDIS archiv^, but no significant change 
was found in fitted position. Given that the different TEC maps make use of many of the 
same GPS stations, this is unsurprising. This problem is exacerbated at southern declinations 
due to the low density of GPS receivers at southern latitudes. Additionally, any errors in 
the TEC maps would have an impact ~ 40% greater for this first epoch, due to its lower 
observing frequency of 1400 MHz, compared with the 1650 MHz center frequency used for 
all subsequent observations. 

Thus, due to the probability of residual ionospheric errors for this epoch, and also the 
potential for frequency-dependent calibrator source structure, the first epoch (MJD 53870, 



^ftp://cddis. gsfc.nasa.gov/pub/gps/products/ionex/ 
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the only 1400 MHz epoch) was dropped from all further analysis. The fit to the remaining 
7 epochs, with ionospheric corrections re-enabled, is shown in Figure [61 While a realistic fit 
is now obtained, the measurement of parallax is still not significant. 



6.3. Data weights 

Initial pulsar imaging and position fitting used visibilities weighted according to the best 
estimate of instantaneous baseline sensitivity. Whilst this is theoretically optimal for data 
which consists only of signal S and additive random thermal errors i?therm, it fails to account 
for the presence of multiplicative systematic errors E^ys = e*''^"^"''*^ caused by residual calibra- 
tion errors. Typically, these systematic errors are dominated by atmospheric and ionospheric 
gradients, although other contributions include antenna and calibrator so urce posit i on er- 
rors, time-variation of antenna bandpasses, and instrumental phase jitter. iFomalontI (120051 ) 



presents a theoretical review of phase referencing errors, while iPradel et al. 
a simulation-based approach. 



J2006h 



presents 



If 0sys(O had zero-mean and was ergodic, its effect would be indistinguishable from 
thermal noise and could be easily estimated, allowing the visibility weights to be corrected. 
For antenna/source position errors and large-scale atmospheric/ionospheric structure, how- 
ever, the residual errors are correlated over long times, causing systematic shifts in the fitted 
position for the target. 

When normal sensitivity-based weighting is employed in the presence of substantial and 
persistent systematic phase errors, the systematic noise on the most sensitive baseline will be 
absorbed into the fit, at the cost of a poorer fit to the less sensitive baselines. The magnitude 
of the induced error will be dependent on the ratio of systematic to thermal errors, and the 
discrepancy in sensitivity between the most and least sensitive baselines in the array. For the 
LBA, the most sensitive baseline (Parkes-DSS43: system equivalent flux density 30 Jy) is 
roughly 13 times more sensitive than the least sensitive (Hobart-Mopra: 380 Jy). Thus, the 
LBA is particularly susceptible to the influence of systematic errors, due to the pronounced 
variation in baseline sensitivities. 

The systematic errors can be crudely estimated (in a model-dependent fashion) by 
performing phase-only self-calibration on the target pulsar over a sufficiently long timescale 
to obtain sufficient SNR, and comparing the magnitude of the corrections to those expected 
from thermal noise alone. While this approach probes systematic errors over a shorter time 
period than those which would dominate for inter-epoch errors (tens of minutes, rather than 
hours to days), it is illustrative of the presence of systematic errors overall. Such corrections 
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are shown in Figure [7] for the ATCA station using a three-minute solution interval during the 
observation on MJD 54127. The corrections are clearly correlated over timescales of tens of 
minutes, and the rms deviation of 4.4 ° is an order of magnitude greater than the estimated 
thermal phase rms of 0.4° (calculated in this high-SNR limit as station sensitivity divided 
divided by target flux, scaled by pulse filtering gain and converted from radians to degrees). 
Thus, for this observation, systematic errors 3> thermal errors and weighting visibilities by 
sensitivity actually degrades the quality of the position fit. 

If the average Egjs could be accurately estimated for each baseline over the duration of an 
experiment, the baseline visibility weights could be adjusted by assuming that i^sys is time- 
invariant. Even more desirable would be the estimation of -Egys as a function of time, allowing 
a time-variable adjustment of the visibility weights. Given present instrumentation, there 
is no way to reliably estimate systematic error (time dependent or independent) in a model 
independent fashion. In the limit where i?sys ^ ^'therm, however, the visibility weight for each 
baseline (regardless of sensitivity) will be dominated by the systematic error contribution, 
resulting in approximately equal weights for all visibilities. Accordingly, visibility weights 
for all baselines were reset to an equal, constant value and data reduction repeated. The 
results are discussed below in §6.4[ 



6.4. Final results 



The revised fit obtained using equally weighted visibility data is described in Table [5] 
and plotted in Figure [HI Through comparison of Table [5] with Table [3], it can be seen that the 
average fit error for a single epoch has increased by 20%, but the inter-epoch systematic error 
has decreased by 95%. Thus, while using equally weighted data incurs a small sensitivity 
penalty, it benefits significantly through the reduced susceptibility to systematic errors. 

The use of natural weighting, as opposed to uniform weighting, was investigated but 
found to produce inferior results. Fitted parameters remained relatively constant but errors 
on the fitted parameters increased by ~ 50%. This is unsurprising, since the use of natural 
weighting promotes a larger beamsize, since more visibility points are concentrated at small 
uv distances. 

The best-fit distance of 26001*1450 consistent with the DM-based distance pre diction 



from the NE2001 Galactic electron dis tribution model (2350 p c JCordes fc Lazidl2002l ). which 
differed considerably from the earlier iTaylor fc Corded (119931 ) distance estimate of 1580 pc. 
Whilst DM-based distance predictions are often assumed to be accurate to ~20%, errors 
up to a factor of several have been observed for individual objects (e.g. PSR B0656+14; 
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Brisken et al.ll2003l ). The VLBI distance is also con sistent with the lower d istance estimate 
of 2.0 ± 0.5 kpc made using HI line absorptions by iKoribalski et al.l (119951 ). The measured 
values of proper motion = 1.52±0.14 mas y r ~^, to = 13.15±0.05 mas yr~^) are consistent 
with the VLA observations of iFomalont et al.l (119971 ). who measured = 1 ± 6 mas yr~^, 
/i5 = 14 ± 11 mas yr~^. Neglecting acceleration from the Galactic potential, the kinematic 
age of the pulsar can be estimated from its Galactic latitude b = 6.3667° and proper motion 
perpendicular to the Galactic plane = 8.93 mas yr~^ as 2.57 ± 0.51 Myr, assuming a 
birth location within 100 pc of the plane. Under the standard assumption of a braking 
i ndex of 3, the obse rved period P = 257 ms and period derivative P = 1.01916 x 10^^^ 
( jSiegman et al.lll993l ) imply a birth p eriod Pn betweeii 35 an d 151 ms - longer than is often 
assumed for normal pulsars (see e.g. iMigliazzo et al.l 120021). but similar to the calculated 
value of Fo = 139.6 ms for PSR ,70538+2817 (Ixramer et al]l2003l ). 



With an accurate proper motion now calculated, the position angle of the proper motion 
for PSR J1559-4438 can be compared to the position angle of the emissi on polarisation. 



which tests the alignment of the rotation and velocity vector, as described by I Johnston et al. 
( 120071 ). If the pulsar emits predominantly parallel or perpendicular to the magnetic field 



lines, then the angle be tween the veloc i ty and polarisation position angles is expected to be 
0° or 90° respectively. iJohnston et al.l (120071 ) found plausible alignment in 7/14 pulsars of 
similar ages to PSR J1559-4438. From Table [5l it is easy to calculate the velocity position 
angle as PA„ = 6.6 ± 0.6°. The polarisation position angle for PSR J1559-4438 given in 



Johnston et al.l (|2007l ) is 71 ± 3, and so in this instance there is no strong case for alignment 
of the velocity and rotation axes. 

One inconsistency with previously published data on PSR J1559-4438 is the transverse 
velocity, which at 164 km s^^ (95% confidence upp er limit of 287 km s~M is inconsistent 
with the 400 km estimated from scintillation by iJohnston et al.l (119981 ). However, this 
scintillation estimate assumes that the scattering material resides in a thin screen at the 
midpoint between the pulsar and the solar system, and neglects any motion of the scattering 
screen itself. As the true distribution of the scattering material along the line of sight to the 
pulsar is not known, the most reasonable interpretation of our results is that the scattering 
screen for PSR J1559-4438 resides considerably closer to the pulsar than to the solar system. 



7. Optimal data weighting 

From the results shown in §6.31 and §6.41 it is clear that for PSR J1559-4438 our astro- 
metric error budget is dominated by systematic errors, and that the use of equally weighted 
visibility points is optimal. However, Table [T] shows that this may not be the case for other 
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pulsars in our target sample, as some are orders of magnitude fainter than PSR J1559-4438. 
Accordingly, we have investigated the conditions under which sensitivity-weighted visibil- 
ities give superior results to equally-weighted visibilities. This was carried out by adding 
simulated thermal noise of varying RMS to the existing dataset. 

Three "noisy" datasets Da, Db, and Dc were constructed by adding gaussian-distributed 
noise to the real and imaginary visibility components of the original observations. Since the 
theoretical single epoch SNR for sensitivity-weighted data should be ~ 800 (a factor of 10 
greater than the typically obtained SNR), the RMS of the added noise in the three datasets 
was set to predicted baseline sensitivity scaled by a factor of 20, 40, and 80, which should 
allow a maximum single-epoch SNR of 40, 20 and 10 respectively. The results of fitting 
the modified datasets (using the inclusive fit approach only) with and without the use of 
sensitivity weighting are presented in Tables Eland [7] respectively. 

Tables [6] and [7] show that while the equally weighted dataset performs better for Da, 
when the average epoch SNR is still high, its performance rapidly deteriorates as the average 
epoch SNR decreases. In Dc, the pulsar was not detected in several epochs using equally- 
weighted data. The reduction in performance is less marked for the weighted datasets, 
although they are clearly still affected by systematic errors. However, if the pulsar was closer 
and the parallax larger, these systematics would be less dominant, and weighted datasets 
would allow measurement of a parallax when equally weighted datasets may be overwhelmed 
by thermal errors, to the point of not detecting the pulsar in a single epoch. 

As noted in §6.31 the use of weighted visibility points would always be optimal if the 
weights could include an estimate of the systematic error contribution to that visibility. 
In the absence of such an estimate, we propose that for the LB A with typical observing 
conditions and calibrator throws, the transition region from systematic to thermal error 
dominated astrometry occurs when the single-epoch detection SNR falls to approximately 
10 for equally-weighted visibilities. This is shown by the similarity of result quality for Db, 
where the average epoch SNR was approximately 10 for the equally-weighted visibilities. 
Alternatively, both weighting regimes can be used and average total single-epoch error (for- 
mal + systematic) can be compared to estimate the optimal weighting scheme. Again, this 
is borne out in the simulated datasets, where a transition from systematic errors dominating 
to epoch fit errors dominating is seen as more noise is added. 
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8. Estimated contributions to systematic error 



The major contributions to systematic error in VLBI astrometry include geometric 
model errors (station/source position, Earth Orientation Parameters), residual ionospheric 
and tropospheric errors, variable phase reference source structure, and image wander due to 
refractive scintillation. Of these, at 1650 MHz residual ionospheric erro rs would be expecte d 
to dominate, despite the a priori ionospheric calibration employed (e.g. iBrisken et al.ll2002l ). 
However, these should be largely uncorrelated from epoch to epoch, and hence the addition 
of more observations can be expected to continue to improve the fit. Residual tropospheric 
errors should also be uncorrelated with epoch, and considerably smaller. 

Geometric model errors cause relative astrometric errors which increase with calibrator 
throw. Earth Orientation Parameters (EOPs) are well d etermined by geod etic observations 
and make minimal contributions to astrometric errors (jPradel et al.ll2006l ). Similarly, the 
mean position of well-determined calibrators makes a minimal contribution. Some LBA 
stations, however, have position uncertainties of several cm, which could make a several 
hundred /las contribution to systematic error. The magnitude of the error depends on 
the source declination and the calibrator-target separation. Given similar uv coverage at 
different epochs, however, the offset will be largely constant with time and is absorbed into 
the reference position of the target. Future planned geodetic observations will continue to 
improve LBA station positions, and reduce this systematic contribution. Additionally, as 
noted in §4.2[ small station position errors can be corrected post-correlation, which offers 
the potential to further improve previous position fits. 

Refractive image wander is caused by large-scale fluctuations in the ISM, and can be 
estimated ba sed on the strength of the pulsar scattering and the scattering disk size (e.g. 
Rickettlll990l ). For strong scattering, which includes PSR J1559-4438, the image wander is 
less than the scatterin g disk size, if Kolmogorov turbulence is assumed for the scattering 
material (jRickettlllQQOl ). Thus, since Table [2] shows that scattering disk of PSR J1559-4438 
is estimated to be only 133 /ias at our observing frequency, the maximum refractive image 
wander is <C 100 /ias, and can be discounted as a source of systematic error. Table [2] shows 
that refractive scintillation is unlikely to be significant for any of our currently targeted 
pulsars. 

The variability of calibrator structure with time depends on the source chosen, but all 
compact extragalac tic radio source s are expected to show some variability, with typical rms 
values of 100 //as (IFomalontl 120051 ). Over short time periods, this image wander may be 
correlated from epoch to epoch and absorbed into proper motion fits, but over long times 
(which could be longer than as astrometric observing program), the mean apparent position 
will be constant. 
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9. Conclusions 

A pipeline for the reduction of LBA astrometric data has been developed in Parsel- 
Tongue and verified by calculating the parallax (vr = 0.384 ± 0.081 mas) and proper motion 
{Ha = 1-52 ± 0.14 mas yi-\ /xg = 13.15 ± 0.05 mas yr"!) of PSR J1559-4438. The cal- 
culated values are consistent with the DM distance estimate and earlier HI absorption and 
proper motion studies. Full account has been made of the impact of residual systematic 
errors on the quality of the astrometric fit. The optimal weighting scheme in the presence of 
systematic errors and varying thermal errors has been investigated, resulting in the guide- 
line that superior astrometric quality can be obtained for the LBA for typical astrometric 
observations by using equally weighted, as opposed to sensitivity weighted, visibilities if the 
target can be detected with S/N > 10. The completion of this parallax program will result 
in the publication of 5 more Southern Hemisphere pulsar parallaxes, which will quadruple 
the number of Southern Hemisphere pulsars with parallaxes determined directly from VLBI 
astrometry. 
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Table 1. Target pulsars 



Pulsar 

name 


DM distance 

(per 


1600 MHz 

flux (mjy) 


Estimated pulsar 

gating gain 


Equivalent gated 
1600 MHz flux (m.Ty) 


Reference 

source 


Calibrator / target 
separation (deg) 


J0108-1431 


130 


0.6 ^ 


5.1 


2.0 


JOlll-1317 


1.5 


J0437-4716 


140 


140 


6.25 


876 


J0439-4522 


1.9 


J0630-2834 


2160 


23 


3.6 


81 


J0628-2805 


0.7 


J0737-3039A/B 


570 


1.6 


2.6 


4 


J0738-3025 


0.4 


J1669-4438 


1600 


40 


3.6 


144 


J1604-4441 


0.9 


J2048-1616 


640 


13 


3.6 


47 


J2047-1639 


0.5 


J2144-3933 


180 


0.8^' 


10 


8 


J2141-3729 


2.1 


J2145-0750 


500 


8>» 


4.3 


34 


J2142-0437 


3.3 



^Taken from Taylor & Cordes (1993) 

^Pulsar suffers heavily from long timescale scintillation, so individual epochs vary considerably from the average value shown. 



Table 2. Observed pulsar scintillation parameters and estimated scattering disk sizes 



Pulsar 
name 


Scintillation 
time (s) 


Scintillation 
bandwidth (MHz) 


Observing 
frequency (MHz) 


Reference Predicted scattering disk^ 
size (mas) at 1650 MHz 


J0108-1431 


*b 


* 


* 


Johnston, Nicastro &: Koribalski (1998) 


0.003 


.10437-4715 


600 


17.4 


660 


Johnston, Nicastro &c Koribalski (1998) 


0.052 


J0630-2834 


466 


0.2 


327 


Bhat, Rao & Gupta (1999) 


0.022 


J0737-3039A 


> 70'= 


1.4=' 


1400 


Coles et al. (2005) 


0.406 


J1559-4438 


77 


0.16 


660 


Johnston, Nicastro & Koribalski (1998) 


0.133 


J2048-1616 


138 


0.54 


327 


Bhat, Rao & Gupta (1999) 


0.023 


J2144-3933 


1500 


2.9 


660 


Johnston, Nicastro & Koribalski (1998) 


0.126 


J2145-0750 


1510 


1.48 


436 


Johnston, Nicastro & Koribalski (1998) 


0.046 


^Calculated from decorrelation bandwidth where available, and taken from NE2001 model where scintillation measurements 
^Scintillation parameters have not been measured, but are believed to be extremely large, as expected for a very nearby pul: 


are unavailable 
sar 



'^Varies considerable with orbital velocity — values quoted for highest transverse velocity 



Coles, private communication 



Table 3. Initial results (sensitivity-weighted visibilities) for PSR J1559-4438 



Parameter 


Naive fit 


Bootstrap fit 


Inclusive fit 


RA (J2000) 


15:59:41.526092 ± 0.000005 


15:69:36.626093 ± 0.000031 


16:69:41.626077 ± 0.000026 


Dec (J2000) 


-44:38:45.902028 ± 0.000017 


-44:38:46.902034 ± 0.000062 


-44:38:46.901989 ± 0.000102 


lice (mas yr~^) 


2.64 ± 0.07 


2.65 ± 0.44 


3.08 ± 0.46 


(mas yr~l) 


13.36 ± 0.02 


13.36 ± 0.06 


13.29 ± 0.17 


7r (mas) 


-0.066 ± 0.038 


-0.054 ± 0.193 


-0.083 ± 0.245 


Nominal distance (pc) 


-15100 


-18500 


-11900 


Nominal vt (km s~^) 


-974 


-1190 


-770 


Reduced chi— squared 


8.8 




1.0 


Average epoch mean fit error(mas) 






0.198 


Average intra— epoch systematic error (mas) 






0.097 


Average inter— epoch systematic error (mas) 






1.070 
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Table 4. Average position shift due to ionospheric corrections for PSR J1559-4438 



Epoch MJD 


RA shift (mas) 


Dec shift (mas) 


Sun separation (deg) 


53870 


-5.36 


0.27 


153.5 


53970 


-1.34 


0.15 


96.4 


54057 


-6.61 


0.39 


25.8 


54127 


-3.15 


-0.07 


63.0 


54182 


-2.66 


0.35 


113.6 


54307 


-2.61 


0.34 


120.1 


54413 


-3.68 


0.20 


30.8 


54500 


-2.89 


0.05 


69.8 



Table 5. Final results (equally weighted visibilities) for PSR J1559-4438 



Parameter 


Naive fit 


Bootstrap fit 


Inclusive fit 


RA (J2000) 


15:59:41.526121 ± 0.000006 


15:59:35.526121 ± 0.000008 


15:59:41.526126 ± 0.000008 


Dec (J2000) 


-44:38:45.901849 ± 0.000020 


-44:38:45.901859 ± 0.000072 


-44:38:45.901778 ± 0.000035 


^ict (mas yr~ ^ ) 


1.62 ± 0.08 


1.60 ± 0.16 


1.52 ± 0.14 


(mas yr~l) 


13.20 ± 0.02 


13.21 ± 0.08 


13.15 ± 0.05 


7r (mas) 


0.280 ± 0.048 


0.291 ± 0.111 


0.384 ± 0.081 


Nominal distance (pc) 


3570 


3440 


2600 


Nominal vt (km s~^) 


225 


217 


164 


Reduced chi— squared 


3.6 




1.0 


Average epoch mean lit error(mas) 






0.242 


Average intra— epoch systematic error (mas) 






0.259 


Average inter— epoch systematic error (mas) 






0.056 



Table 6. Noise-added fits for PSR J 1559-4438, sensitivity-weighted visibilities 



Parameter 


Da 






Db 




Dc 






RA (J2000) 


15:59:41.526163 


± 


0.000027 


16:69:35.526086 ± 0.000025 


16:69:41 


,526025 


± 


0.000041 


Dec (J2000) 


-44:38:45.902054 


± 


0.000154 


-44:38:45.902131 ± 0.000163 


-44:38:45 


.901760 


± 


0.000258 


/Xq, (mas yr""*^) 


2.03 


± 


0.46 


1.44 ± 0.59 




2.83 


± 


0.79 


(mas yr~^) 


13.39 


± 


0.24 


13.40 ± 0.21 




13.00 


± 


0.44 


TT (mas) 


0.109 


± 


0.278 


0.270 ± 0.373 




0.621 


± 


0.579 


Nominal distance (pc) 


9200 






3700 




1610 






Nominal vt (km s""'") 


690 






237 




102 






Reduced chi— squared 


1.0 






1.0 




0.5 






Average epoch mean fit error(mas) 


0.441 






0.904 




1.308 






Average intra— epoch systematic error (mas) 


0.604 






1.412 




0.863 






Average inter- epoch systematic error (mas) 


0.769 






0.315 




0.0 






Average single— epoch SNR 


33 






17 




9 
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Table 7. Noise-added fits for PSR J1559-4438, equally weighted visibilities 



Parameter 


Da 


Db 




Dc 






RA (J2000) 


15:59:41.526142 ± 0.000026 


16:69:35.526115 ± 0.000042 


16:69:41 


.526268 


± 


0.000084 


Dec (J2000) 


-44:38:45.901745 ± 0.000076 


-44:38:45.902220 ± 0.000234 


-44:38:45 


.901798 


± 


0.000420 


fj-cr. (nias yr""*") 


0.91 ± 0.38 


0.22 ± 0.99 




-0.04 


± 


1.20 


fj.S (mas yr~^) 


13.18 ± 0.10 


13.56 ± 0.32 




12.65 


± 


0.65 


TT (mas) 


0.544 ± 0.235 


0.760 ± 0.597 




2.092 


± 


0.787 


Nominal distance (pc) 


1840 


1320 




480 






Nominal vt (km s~ ) 


116 


85 




29 






Reduced chi— squared 


0.8 


1.0 




0.9 






Average epoch mean fit crror(mas) 


0.618 


0.937 




1.853 






Average intra— epoch systematic error (mas) 


0.791 


1.433 




1.531 






Average inter— epoch systematic error (mas) 


0.000 


0.380 




0.000 






Average single— epoch SNR 


20 


10 




5 









U (10* X) 



U (10° X) 



Fig. 1. — Typical uv coverage at 1650 MHz (no Ceduna; left panel) and at 8400 MHz (with 
Ceduna; right). The target sources are PSR J 1559-4438 and PSR J0437-4715 respectively. 
Without Ceduna, the uv coverage is heavily biased North-South, and wide hour-angle cover- 
age is necessary to gain acceptable uv coverage. The uniformly weighted beam size is 40x13 
mas at 1650 MHz, and 3.2x2.6 mas at 8400 MHz. 




Fig. 2. — The effects of diffractive scintillation for pulsar J1559-4438. Panel a) shows the 
uncorrected visibility amplitude on baselines including the Tidbinbilla antenna from one 
experiment - the scintillation timescale of several minutes is apparent. Panel b) shows 
the same visibility amplitudes after correcting for scintillation. Panel c) shows the image 
residuals of the uncorrected dataset - contours are 1,2,4 and 8 mJy/beam. Panel d) shows 
image residuals after correcting for scintillation, with contours at 1,2 and 4 mJy/beam - the 
improvement in image quality is obvious. 
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Fig. 3. — Motion of the pulsar in right ascension and dechnation, with measured positions 
overlaid on the best fit. Sensitivity-weighted visibilities were used. The motion of the 
pulsar is positive in right ascension and declination. The first epoch (lower left) is clearly 
inconsistent. 
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Fig. 4. — Variation of fitted position shifts due to ionospheric correction for PSR J1559-4438 
vs angular separation between the pulsar and the Sun. The single 1400 MHz epoch has an 
unusually large correction despite the large angular separation between the pulsar and the 
Sun. 
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Fig. 5. — Motion of the pulsar in right ascension and dechnation, with measured posi- 
tions overlaid on the best fit, when ionospheric corrections have been removed. Sensitivity- 
weighted visibilities were used. The first epoch (lower left) is now more consistent, but the 
third epoch (during which the pulsar-Sun angular separation was only 26 °) is now inconsis- 
tent. 
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Fig. 6. — Motion of the pulsar in right ascension and dechnation, with measured positions 
overlaid on the best fit, with ionospheric corrections reinstated but the first epoch dropped. 
Sensitivity-weighted visibilities were used. The fit is improved considerably. 



-29- 



150- 



100 



V) 



c 
o 

■l-l 

u 

o 
u 




■100 



-150 



0.75 0.80 0.85 0.90 0.95 1.00 

Time (Days since MJD 54127) 



1.05 1.10 



Fig. 7. — Self-calibration corrections, using a three-minute timescale, for the ATCA station 
on J1559-4438. Clear systematic deviations are seen from the zero-mean distribution ex- 
pected from purely thermal noise. The rms of the corrections exceeds those expected due to 
thermal noise by an order of magnitude. 
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Fig. 8. — Motion of the pulsar in right ascension and dechnation, with measured positions 
overlaid on the best fit, with ionospheric corrections reinstated but the first epoch dropped. 
Equally weighted visibilities were used, mitigating systematic errors and allowing for the 
first time a significant measurement of the parallax for this pulsar. 



